An SDOF lightly damped viscoelastic system with right unilateral nonzero offset barrier impacts behavior under random excitations is investigated, which is schematically shown in Fig. 1. The equation of the system is of the form 1a 1b

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic model of the viscoelastic oscillator with right unilateral nonzero offset barrier impacts.

In Eq. (1a), ε is a small positive parameter, and the dot indicates differentiation with respect to time; represents lightly linear or nonlinear damping force; r(x) denotes the nonlinear or linear restoring force as an odd function of x; η_m (t) represents Gaussian white noises in the sense of Stratonovich, and possesses correlation functions E[η_n(t)η_k (t + τ₀)] = 2D_nk δ (τ₀); and (m ∈ N⁺) describes amplitudes of weakly random perturbations. Equation (1b) is the impact condition for the inelastic instantaneous interaction, where ρ in [0,1] is a restitution factor. and represent the velocity of the system before and after the instant of impacts; S₀ denotes the position of the nonzero rigid barrier.

The constitutive model of viscoelastic damping is provided as^[11] 2 and the relaxation function q(t) is expressed as follows based on Refs. [9] and [11]: 3 Each λ_i in Eq. (3) is the relaxation time of the viscoelastic components, and β_i specifies elastic modulus being either positive or negative. Both parameters are determined by the concrete problems.

Taking into account that the viscoelastic force contributes to both the stiffness and the damping,^[16] system equations (1a) and (1b) can be written as 4a 4b Here, the consideration of the viscoelastic force has been incorporated in and r₁ (x), respectively. The potential energy function U(x) and the total energy E of the system are 5 6 The period of the corresponding free motion is given by 7 where A is the amplitude for a given E, and the positive root is calculated from E = U(A).

The motion equations of the system (4a) and (4b) can be described through variable substitution. Utilize the total energy E(t) and the total phase angle ϑ(t) according to the transformation 8 and ϑ(t) is defined as 9 ϕ (t) denotes the residue phase angle, and ω(t) is the instantaneous frequency. The average frequency of ω(t) can be calculated by the quasi-period T 10 Denote each component in the viscoelastic force exhibited in Eqs. (2) and (3) as 11 Applying the transformation of Eq. (8) with , and neglecting the decaying transient part, we can get 12 Thus, functions and r₁ (x) in Eq. (4a) are derived by Eq. (12) as 13 Besides, combining Eqs. (7) and (10), the average frequency is determined 14 Thus, viscoelastic system equations (1a) and (1b) with unilateral nonzero offset barrier impacts are transformed to the system Eqs. (4a) and (4b) without viscoelastic damping form.

Introducing , equation (4a) can be rewritten as the following Itô stochastic differential equations: 15 where Γ(x,y) and Σ(x,y) denote drift and diffusion representations, and W(t) is a unit Weiner process.

The potential energy function and total energy of the system have been reported in the previous Eqs. (5) and (6), respectively. When the parameter ε is zero, equations (15) and (4b) constitute a free VI system. Suppose functions r₁ (x) and U(x) are such that the free VI system associated with Eq. (15) without impact has a gens periodic solutions surrounding the origin on phase plane (x,y). The period is the same as earlier 16 Here, a denotes the maximum displacement derived by E = U(a). As E > U(S₀), the system will impact the nonzero offset barrier. The system state before impact is , and the system energy can be written as 17 When the impact occurs, the system state after impact with the nonzero offset barrier turns into , and the system energy becomes 18 The loss of the system energy during each impact is calculated as 19 The system still moves as a free oscillator between the first impact (t = t₁) and the second impact (t = t₂), and the time period Δt of the process between both consecutive impacts can be evaluated by 20 The time period can be seen as the period of free oscillation after occurrence of impacts, which relies on the system energy; thus, denoted as T(E) = Δt represents the quasi-period. It is implied from Eq. (19) that the energy loss ΔE is related to the restitution factor ρ, the energy before nonzero rigid barrier impacts, the nonzero position S₀, and the viscoelastic force due to its equivalent stiffness part.

Relying on the Eqs. (6) and (15), the Itô equations for displacement x and system energy E can be derived by employing the Itô differential rule. When the damping, the stochastic excitation, and the energy loss are small, it can be found that energy E(t) is a slowly varying process, and the displacement x(t) is a rapidly varying process. Thus, the method of SAEE^[45] is absorbed here. Further, the averaged Itô equation for E(t) is given as 21 where Γ(E) and Σ(E) are averaged drift and diffusion representations, and B(t) is the standard Brownian motion process.

For the situation of E < U(S₀), the system keeps the vibration without impact. The drift and diffusion representations of the viscoelastic system with nonzero offset barrier impacts remain as 22a 22b 22c As E exceeds U(S₀), the system possesses the vibration with the impact. The energy loss during impact needs to be counted; thus, the averaged drift and diffusion representation becomes 23a 23b 23c where A = U⁻¹(E) denotes the amplitude of vibration, and ΔE[T(E)]⁻¹ describes the energy loss during each impact.

Further, the averaged Fokker–Plank–Kolmogorov (FPK) equation associated with Itô Eq. (21) is expressed as 24 where p = p(E, t|E₀) represents the transition probability density of E dependent on the initial condition 25 and the boundary condition is finite for the situation of E = 0, while 26 Thus, the stationary solution of Eq. (24) can be provided by 27 where Z is the normalization constant.

Next, the joint probability density of the displacement and the velocity can be calculated 28 Here, p(x|E) is the conditional probability density function, which can be given by 29 The normalization constant Z′ is computed by 30 Substituting Eqs. (29) and (30) into Eq. (28) yields 31 Then, the marginal probability density functions p(x) and p(y) can be solved by 32a 32b

Here, we adopt the nonlinear viscoelastic oscillator with right nonzero rigid barrier impacts to illustrate the analytical procedure 33 The damping force and viscoelastic force V′ here are assumed to be of order ε, and η(t) are stationary broad-band random processes of order ε^1/2, i.e., the Gaussian white noise with intensity 2D′. σ denotes the nature frequency. ρ and S₀ are described the same as the previous interpretation. By approximating V′ based on Eqs. (12) and (13) to the case, the equivalent system of Eq. (33) has the form 34 where and .

According to Eq. (14), specifically for the present case, the average frequency can be given by 35 Relying on the method of SAEE, the drift and diffusion representations of the averaged Itô equation for energy E can be obtained in detail.

For the situation of E < U(S₀), combining with Eqs. (22a)–(22c), let , and then yield 36

As E exceeds U(S₀), combining with Eqs. (23a)–(23c), we can get 37

Next, by solving the averaged FPK Eq. (24), we obtain the PDFs of the system energy E, and further get the probability density functions p(x,y), p(x), and p(y) based on Eqs. (31), (32a), and (32b). And then consider the case of , system equation (33) represents a nonlinearly damped and stochastically excited Van der Pol viscoelastic VI system.

The analytical stationary PDFs will be exhibited and compared with the directly numerical simulation results. All the directly numerical simulations computed by employing the Runge–Kutta (RK) algorithm and Monte Carlo method are presented based on the original system Eq. (33), which provide the validation as well as rationality of the analytic results. Solutions are initially observed with the parameters: d₀ = −0.01, d₂ = −0.005, σ = 1.0, β₁ = −0.2, λ₁ = 1.1, ρ = 0.9, S₀ = 0.995, and D′ = 0.0236, unless otherwise stated. Then, the effects on the system responses for the intensity D′ of stochastic process η(t), two viscoelastic parameters β₁ and λ₁, the position S₀ of the right nonzero rigid barrier, and the restitution factor are described in detail in Figs. 2–5. The increase of noise intensity D′, i.e., gaining stronger external force leads to lower peaks of stationary PDFs of the displacement x and the velocity y, and smaller values for the stationary PDFs of total energy E at small energy, i.e., make ones more flat in Figs. 2(a)–2(c) and the oscillator stays at equilibrium with a smaller possibility. Besides, figure 2(a) shows that the increase of D′ can contribute to the system with high energy, large displacement, and velocity, i.e., obtaining larger probability of the system. It can be seen that the approximate analytical results are relatively close to the directly numerical ones.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Stationary PDFs of (a) total energy E, (b) displacement x, and (c) velocity y, solved both analytically and numerically with different noise intensities. Solid lines represent analytical results, while circles, squares, and triangles represent Monte Carlo simulation results.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Stationary PDFs of (a) total energy E, (b) displacement x, and (c) velocity y, solved both analytically and numerically with different viscoelastic parameters. Solid lines represent analytical results, while circles, squares, and triangles represent Monte Carlo simulation results.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Stationary PDFs of (a) total energy E, (b) displacement x, and (c) velocity y, solved both analytically and numerically with different positions of the nonzero rigid barrier. Solid lines represent analytical results, while circles, squares, and triangles represent Monte Carlo simulation results.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Stationary PDFs of total energy (a) E, (b) displacement x, and (c) velocity y, solved both analytically and numerically with different restitution factor values. Solid lines represent analytical results, while circles and squares represent Monte Carlo simulation results.

The elastic modulus β₁ and relaxation time λ₁ show a difference in responses of the system according to Fig. 3. β₁ keeps negative in Fig. 3, and the change of stationary PDFs caused by the variation of |β₁| is stronger than that caused by the variation of λ₁. The increase in |β₁| is equivalent to an increase in the damping effect, which can contribute to higher peaks of stationary PDFs of displacement and velocity, and larger values for the stationary PDFs of total energy E at small energy, indicating a larger probability that the system will stay close to the equilibrium state. The two solutions here also agree well.

Examination of the position S₀ of the right nonzero rigid barrier and the restitution factor ρ on the system response are illustrated in Figs. 4 and 5, respectively. Curves in Fig. 4(a) show unsmoothness at E = U(S₀). When the oscillator is further away from the obstacle position in Figs. 4(b) and 4(c), which is equivalent to reducing the damping effect, it leads to lower peaks of stationary PDFs of displacement and velocity, and smaller values for the stationary PDFs of total energy E at small energy, indicating a smaller probability that the system will stay close to the equilibrium state. Additionally, when S₀ reaches the amplitude of the oscillator motion, i.e., the maximum displacement which is approximately 1.995 here, the corresponding response probability is close to 0, indicating that if S₀ crosses the threshold, the impacts will not occur. In Figs. 5(a)–5(c), the variation of restitution factor ρ from 0.45 to 0.94 does not cause a great change in the response probability of the system. It can be found that the analytical results and MC simulation results agree well in these figures, which further indicates that the analytical procedure is reliable and effective via the related comparisons.

According to Section 4.1, satisfactory agreement is shown between the results of analytical procedure and MC simulation solution obtained from the original stochastic VI system. Based upon the above conclusion, a further analysis for joint stationary PDFs of the displacement x and the velocity y is observed here, and stochastic P-bifurcation occurs with the variation of certain parameters. Thus, we perform the exploration for stochastic P-bifurcation in this subsection. Take ρ = 0.95 and S₀ = 0.985, and the rest are kept the same as the initial values.

The position S₀ of the nonzero rigid barrier varies from 0.985 in Figs. 6(a), 6(b), and 6(d) to 1.985 in Fig. 6(c), and the figures exhibit the decreasing joint PDFs and the corresponding analytical and numerical contour lines for different elastic moduli β₁ being negative, which is the response to the result in Fig. 3. The contours of the two solutions further change from Fig. 7(a) to Fig. 7(d) with the large-amplitude oscillations and the small-amplitude oscillations around the origin at the position of S₀ = 0.485 when some parameters are in the bifurcation schemes, and the phenomenon here will be discussed below. Besides, the analytical solutions and MC simulation ones on the contours of joint PDFs agree well here, which responds to the previous conclusion.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The contours for joint PDFs of the displacement x and the velocity y solved both analytically (solid line) and numerically (pecked line) with different elastic modulus values β1: (a) β1 = −0.5, (b) β1 = −0.19, (c) β1 = −0.19, S0 = 1.985, and (d) β1 = −0.06.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The contours for joint PDFs of the displacement x and the velocity y solved both analytically (solid line) and numerically (dotted line) with different bifurcation parameters: (a) β1 = 0.0, (b) λ1 = 0.01, (c) ρ = 1.0, and (d) d0 = −0.05.

Next, figures 8, 10, 12, and 14(a)–14(c) provide joint stationary PDFs of the displacement x and the velocity y and the section graphs of the joint PDFs on the surface x = 0.485 for different damping coefficients d₀, elastic moduli β₁, relaxation times λ₁, and restitution coefficients ρ, respectively. When ρ = 0.95, β₁ = −0.1, and λ₁ = 0.5 are selected, the joint PDFs decreases as |d₀| increases, reflecting the effect of the reduced damping on the small displacement and velocity. In Fig. 8(c), it demonstrates the process of evolving from a unimodal shape of joint PDFs, which shows at the position of the rigid barrier, to a bimodal shape, i.e., the shape of crater in Fig. 8(b). Figures 9(a)–9(d) further provide the corresponding contours. When d₀ = −0.015 in Fig. 9(a), there exists one stochastic attractor, characterized by small-amplitude oscillations around the origin. And when d₀ reaches −0.035, the stochastic attractor is transformed into the characterization of large-amplitude oscillations. When d₀ is further changed to −0.045 and −0.05 in Figs. 9(c) and 9(d), the small-amplitude oscillations around the origin appear again and get stronger, while the existing large-amplitude oscillations is weakened. Clearly, it indicates that the occurrence of stochastic P-bifurcation exists in the process of variations of d₀ based upon the concept of stochastic bifurcation (phenomenological or P-bifurcations are observed when the underlying probabilistic structure of the long-term behavior of the state variables undergoes topological changes with certain parameter variations),^[46–49] and it occurs when d₀ is close to −0.029, and the stability of system is changed here.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. The surface for joint PDFs of the displacement x and the velocity y with the damping coefficient (a) d0 = −0.015 and (b) d0 = −0.05; (c) the sectional views of joint PDFs on the surface x = 0.485 with variations of the damping coefficient d0.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Contour plots for joint PDFs of the displacement x and the velocity y: (a) d0 = −0.015, (b) d0 = −0.035, (c) d0 = −0.045, and (d) d0 = −0.05.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. The surface for joint PDFs of the displacement x and the velocity y with the elastic modulus (a) β1 = 0.0 and (b) β1 = 0.06; (c) the sectional views of joint PDFs on the surface x = 0.485 with variations of the elastic modulus β1.

When d₀ returns to the initial value, a non-negative value of the elastic modulus β₁ of viscoelastic force induces a similar phenomenon. Figures 10(a) and 10(b) illustrate the joint stationary PDFs with the shape of crater for β₁ = 0.0 and β₁ = 0.06, which are different from the effect in Fig. 6. It intuitively shows the evolution process of the joint PDFs from the unimodal shape to the shape of crater caused by varying β₁ from negative to zero and then to positive (equivalent to reducing damping) through the section graphs of joint PDFs when x reaches the position of rigid barrier in Fig. 10(c). Hence, it can be seen that the stochastic P-bifurcation occurs when β₁ approximates to −0.02, which is treated as the watershed of joint PDFs from the singular peak to the shape of crater. Figures 11(a)–11(d) also report the corresponding contours. The variation of the stochastic attractor caused by varying β₁ from negative to zero and then to positive is close to the effect of the increase of |d₀| in Fig. 9.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Contour plots for joint PDFs of the displacement x and the velocity y: (a) β1 = −0.21, (b) β1 = 0.0, (c) β1 = 0.03, and (d) β1 = 0.06.

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. The surface for joint PDFs of the displacement x and the velocity y with the relaxation time (a) λ1 = 0.01 and (b) λ1 = 2.7; (c) the sectional views of joint PDFs on the surface x = 0.485 with variations of the relaxation time λ1.

Further, we examine the relaxation time λ₁ of viscoelastic force from the vicinity of the usual value 1.0 to a wider range, and obtain interesting observations. Figures 12(a) and 12(b) exhibit joint stationary PDFs possessing the shape of crater for λ₁ = 0.01 and λ₁ = 2.7, which are different from the effect in Fig. 3. Figure 12(c) intuitively shows the evolution process of the joint PDFs from the unimodal shape to the shape of crater as λ₁ varies from 0.01 to 0.37 and from 1.59 to 2.7 (equivalent to reducing damping) through the section views of joint PDFs when x reaches the position of rigid barrier. Clearly, it can be seen the stochastic P-bifurcation occurs when λ₁ approximates to 0.13 or 1.87. The corresponding contours depicted in Figs. 13(a)–13(d) show that the transformation on the characterization of the stochastic attractor from small-amplitude oscillations around the origin to large-amplitude oscillations indeed exists in the above two intervals, which has not been covered before.

Fig. 13.

	Figure Option ViewDownloadNew Window
	Fig. 13. Contour plots for joint PDFs of the displacement x and the velocity y: (a) λ1 = 0.01, (b) λ1 = 0.5, (c) λ1 = 1.5, and (d) λ1 = 2.7.

Fig. 14.

	Figure Option ViewDownloadNew Window
	Fig. 14. The surface for joint PDFs of the displacement x and the velocity y with the restitution coefficients (a) ρ = 0.917 and (b) ρ = 1.0; (c) the sectional views of joint PDFs on the surface x = 0.485 with variations of the restitution coefficients ρ.

Then choose β₁ = −0.1 and λ₁ = 0.65, we further consider the effect of the restitution factor ρ limited to the range 0.9–1.0 on the system response. The effect of impact is to increase the damping associated with the system, and qualitatively, decreasing ρ is equivalent to increasing the damping parameter. The stationary joint PDFs are computed for the cases ρ = 0.917 and 1.0, as shown in Figs. 14(a) and 14(b). It decreases for the small displacement and velocity when ρ reaches 1.0, representing the classical elastic impact, and the shape of the crater also appears in the figure. The section views of joint PDFs in Fig. 14(c) further describe the evolution process caused by ρ on the surface of rigid wall and report the watershed quickly, i.e., when ρ approximates to 0.968, the stochastic P-bifurcation occurs, and it becomes double peaks when ρ = 1.0, corresponding to the shape of the crater in Fig. 14(b). Figures 15(a)–15(d) also provide the corresponding contour plots. The characterization of the stochastic attractors develops into large amplitude oscillations at ρ = 0.986 in Fig. 15(b), and both small amplitude oscillations and large amplitude oscillations are visible in Figs. 15(c) and 15(d).

Fig. 15.

	Figure Option ViewDownloadNew Window
	Fig. 15. Contour plots for joint PDFs of the displacement x and the velocity y: (a) ρ = 0.917, (b) ρ = 0.986, (c) ρ = 0.9977, and (d) ρ = 1.0.